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Abstract 


A linear  algorithm  1h  given  for  the  generation  of  covariance 
sequences  for  rational  digital  filters  using  numerator  and  denominator 
coefficients  directly.  There  is  no  need  to  solve  a Lyapunov  equation 
or  to  solve  for  the  residues  of  a spectrum,  as  in  other  methods.  By 
appealing  to  certain  results  from  the  theory  of  inners  we  show  the 
algorithm  provides  a unique  solution  provided  only  that  the  filter 
is  stable. 

Our  results  may  be  used  to  compute  error  variances  due  to  product 
rounding  and  signal  quantization,  and  to  generate  covariance  strings 

used  in  other  studies  involving  second-order  properties  of  digital 
filters. 


I.  Introduction 


There  are  numerous  problems  in  the  study  of  digital  filters  that 
require  solution  of  an  integral  like 


nrr  i huw*-1) 


l,dz 

z 


-0  2n  j 
or,  more  generally,  like 


rk  " 2SJ  * H(z)H(z_1)7*k 


(i) 


(2) 


For  example,  in  the  computation  of  steady-state  output  noise  variance 
o due  to  fixed-point  quantization,  one  needs  o - (A  /12)r0,  where  A 
is  the  quantization  step  size  and  H(z)  is  the  transfer  function  between 
the  quantization  noise  source  and  the  filter  output.  The  same  kind  of 
computation  is  required  for  the  computation  of  output  noise  variance 
due  to  input  signal  quantization.  More  generally,  sequences  like 

K 

{r^}^  are  required  in  the  study  of  model  reduction  procedures  for 
approximating  high-order  filters  with  low-order  ones. 

The  most  obvious  way  to  compute  r^  is  to  perform  the  indicated 
contour  integration  by  evaluating  residues.  Mitra,  et.  al.  [2]  have 
tabulated  the  appropriate  residues  for  generic  terms  arising  in  the 
computation  of  r^.  With  some  tedium  the  results  may  be  generalized 
to  handle  the  calculation  of  r^.  To  apply  the  results  of  [2],  one 
must  solve  for  poles  of  H(z).  This  can  be  a nuisance,  so  one  looks 
for  alternatives. 


Note  r^  may  be  written  [1] 


rk  " h0hk  + c'*k,C  c 


k > 0 


(3) 


where  K is  the  solution  to  a matrix  Lyapunov  equation,  $ is  a state 

CD 

transition  matrix,  {h^l^  is  the  unit  pulse  response  sequence  for  H(z), 


and  c'  - (1  0...0).  The  computational  difficulty  with  (3)  is  that  a 

Ic 

Lyapunov  equation  must  be  solved  and  ♦ determined.  The  latter  deter- 
mination involves  solving  for  the  eigenvalues  of  4>,  performing  N 
matrix  multiplies  to  initialize  the  recursion  to  follow,  or  exploiting 
special  properties  of  the  matrix  ♦ . Another  alternative  is  to  proceed 
as  Jury  does  [3]  and  write  r^  as  the  ratio  of  two  determinants.  This 
alternative  Is  not  particularly  attractive  because  it  does  not  easily 
generalize  to  the  computation  of  r^. 

Our  approach  follows. 

II.  An  Alternative  Method 

Let  H(z)  denote  a stable  autoregressive  moving  average  digital 
filter  of  the  form 


H(z)  - 


M 

I b z' 
m-0  m 


N 

T.  a.z 
. - 1 
1-0 


-1 


5 ao  - 1 


a^,  b^  real  numbers 


(4) 


-1 


N 

I a z ~ - 0 ► |z|  «■'  1 
1-0 


Denote  this  filter  H(z)  : ARMA(N,M).  There  is  no  need  for  M to  be  less 


than  N.  The  covariance  sequence  associated  with  H(z)  is  {r.  } with 

k 

\ - r_k  given  by 

rk  ' zTf  * s<l)zk  T • ,k 


(5) 


S(z)  - H(z)H(z_1) 


Here  the  contour  C lies  within  the  region  of  absolute  convergence  of 
H(z).  The  contour  may  be  chosen  to  be  the  unit  circle  in  which  case 
S(r.  - exp(j2irf>)  is  the  spectrum  (or  magnitude-squared  frequency 
response)  corresponding  to  H(z). 
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ao 

The  unit  pulse  sequence  corresponding  to  H(z)  is  {h.  with 


0 , k < 0 


N 

b.  - Z a h. 
k , n K-n 
n-1 


(6) 


k > 0 


Here  is  assumed  zero  for  k >_  M+l.  The  covariance  sequence  is 
related  to  the  unit  pulse  sequence  as  follows: 


r.  - Z h h ,, 
k . n n+k 
n-0 


r.  - r . 
k -k 


k > 0 


(7) 


Substitute  (6)  into  (7)  to  get 


N 

Z a r.  * d, 
n-0  n k_n  k 


M-k 

i.  - Z hb  . 
k „ n n+k 
n-0 


k > 0 


(8) 


Note  d^  - 0 for  k _>  M+l,  in  which  case  the  {r^}  sequence  behaves  just 
like  a purely  autoregressive  one.  That  is,  the  sequence  {r^l  satisfies 
a linear  homogeneous  difference  equation  for  k > M + 1. 

From  (8)  it  is  clear  that  the  covariance  sequence  may  be  generated 
recursively.  The  trick  is  to  initialize  the  recursion  by  finding  r^ 
for  0 <_  k < N.  Write  out  (8)  for  k - 0,1,... N: 

Ar  - 9 


Va0 

a3+ai 

Va2 


Vao 

a5+al 


®N  aN-l  * 


«5  . . 
Vao 


“N 

“N  ° 
aN  0 . 


ai  ao 


(9) 


4 


r'  - (rQ,  r^...,^)  d'  - (dQ,  dlt...,dN) 

The  matrix  A is  generated  as  follows:  begin  with  the  first  row 
(aQ,  ax,  a2,...,&N);  left-shift  this  row  (n-1)  times  and  add  the  mth 
overflow  to  the  (n,m+l)  term  to  get  the  n^  row.  Thus  the  (i,j)  element 
of  matrix  A is 


IA> 


l.J 


ai-l  for  J “ 1 


Vj-2  + *1-1  tor  J ’ 1 


(10) 


where  a^,  a^t...,an  are  defined  in  (4)  and  a{  ■ 0 for  i < 0 and  i > N. 

For  example,  the  fourth  row  Is  (a^,  a^+a2,  a5+aj»  a6+a0* ’ * ’ ,aN’ 

See  Appendix  I of  |4].  Tlie  right-most  triangular  region  of  the  matrix 

consists  of  zeros. 

M 

The  Method.  The  solution  method  is  to  generate  {h^ig  from  (6),  solve 

_ N M 

(9)  for  r,  and  then  use  (h^),  (a^Jg  ’ an<*  ^bm^0  to  8enerate  an  arbi- 
trarily long  finite-length  version  of  (r^)^,  from  (8).  If  only  (r^}^ 
is  required  (as  in  applications  requiring  only  r^) , then  only  (9) 
must  be  solved  using  standard  techniques  for  solving  linear  equations. 
The  same  procedures  were  outlined  in  [4]  for  generating  covariance 

sequences  for  autoregressive  filters.  For  autoregressive  filters, 

■>  , 

‘*0  " b()  am*  l*k  * 0,  k N 0.  Tlie  matrix  A has  also  arisen  In  Jury  s 

work  15 | . 


Ill,  Comments  on  Unicity 
N M 

The  connection  between  the  (a  }_  and  {b  }.  parameters  of  an  ARMA(N,M) 

n u mu 

00 

filter  H(z)  and  the  corresponding  covariance  sequence  (r^)_(p  is  unique, 
provided  H(t)  is  stable  and  minimum  phase  and  there  are  no  pole-zero 
cancellations.  If  H(z)  is  stable  and  there  exist  no  pole-zero  cancella- 
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tlous,  but  H(z)  la  not  minimum  phase,  then  corresponding  to  H(z)  la  a 
unique  {r^l*^  • but  not  vlce-veraa.  Thus  given  the  feedback  and  feed- 
forward parameters  {a^l^  anc*  f°r  a stabl«  N(z)  one  may  solve  for 

a unique  r from  (9)  and  recursively  get  the  corresponding  unique  covari- 
ance sequence  {r^)”^  . The  following  theorem  is  arcane  and  relevant. 
Theorem:  The  Matrix  A ia  Nonsingular  for  Stable  Filters.  Let  H(z) 

be  the  transfer  function  of  a stable  ARMA(N,M)  filter,  as  defined  in 

N 

(4).  That  is,  the  roots  of  A(z)  » T.  a z lie  strictly  inside  the 

n“0  n 

unit  circle  |z|  - 1.  Then  the  matrix  A of  (9)  is  nonsingular  and 
the  solution  for  r in  (9)  is  unique. 

Proof:  Rotate  the  matrix  A by  a/2  and  interchange  the  order  of  columns 
to  obtain  the  matrix 


ao 

al 

a2 

a3 

aN 

al 

a,+a0 

a3+al 

• • • 

aN-l 

* 

*2 

a3 

a4+a0 

| • • • 

• 

0 

> 

0 

0 

0 ...  o 

a0  _ 

preserve  |il 

- |A| 

. Define 

A.  (z) 

N 

- z 

A(z) 

(11) 


N 

E a z 
n«0  n 


N-n 


N 

1!  i»  z 

n n 
n»() 


n 


(12) 


u - a„ 
n N-n 

The  roots  of  A^z)  lie  inside  |*|  •=  1 by  virtue  of  the  fact  that  the 
roots  of  A(z)  do.  Jury  [5]  shows  that  necessary  and  sufficient  con- 
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dltions  for  the  roots  of  A^z)  to  lie  Inside  |z|  - 1 are 
i)  A1(+l)  > 0 , (-1)°  A1(-l)  > 0 

11)  ■ xn_i  'Yn-1  Poslt*vc  lnnerwise 

In  order  to  relate  the  determinant  of  fl  to  this  stability  criterion, 
define 


"a  a , 

a -» 

* * * 

”o  0 ... 

a" 

n n-1 

n-2 

0 

0 

0 a 

a , 

a. 

an 

a. 

n 

n-1 

1 

• Y ■ 

* n+1 

0 

1 

0 

a 

• 

n 

• 

• 

ao 

0 

a 

a 

n 

0 

n 

- 

m 

• 

(13) 


It  follows  from  Jury  [5]  that 


2I«I  " lXn+l+Yn+l'  " l<-l>nA1(l>Aiei)An-1l  °4) 

If  H(z)  is  stable,  it  follows  that  A^  (±1)  4 0,  4 0,  and 

therefore  that  |ft|  - |a|  4 0.  Q.E.D. 

The  following  second-order  example  is  illustrative. 

Example:  Assume  H(z)  is  second  order.  Then 


This  matrix  is  nonsingular  provided  a^  4 1 and  a^  4 t^+l) . These 
conditions  are  illustrated  in  Figure  1.  Thus  A Is  nonsingular  on  the 
two  dimensional  plane  (a^,  a^),  minus  the  boundary  lines  illustrated. 
Note  the  interior  of  these  lines  is  the  region  of  stability  for  a 
second-order  filter.  Thus  for  a stable  second-order  filter  the  matrix 


7 

A Is  nonsingular.  The  implication  does  not  go  the  other  way.  All 
the  (a^,a2>  pairs  outside  the  illustrated  boundaries  also  give  nonsingu- 
lar A,  but  unstable  H(z).  Solution  of  (9)  for  such  pairs  gives  {rk> 
sequences  that  are  not  covariance  sequences. 

IV.  Remarks 

The  results  presented  here  may  be  used  to  compute  error  variances 
in  finite  word-length  digital  filters  and  to  generate  covariance  strings 
*rk*0  for  U8e  in  model  reduction  procedures  and  the  like.  The  computa- 
tions seem  simpler  than  those  of  [6]  and  generalize  the  results  of  [2]- 
[4],  Software  is  available  from  the  authors  upon  written  request. 
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